Lab Session 10

autoplot(usnetelec)

No transformation required

autoplot(mcopper)

(lambda <- BoxCox.lambda(mcopper))

## [1] 0.1919047

mcopper %>% BoxCox(lambda=lambda) %>% autoplot

autoplot(enplanements)

(lambda <- BoxCox.lambda(enplanements))

## [1] -0.2269461

# I don't like such strong transformations. Will use 0 instead
enplanements %>% BoxCox(lambda=0) %>% autoplot

autoplot(a10)

(lambda <- BoxCox.lambda(a10))

## [1] 0.1313326

a10 %>% BoxCox(lambda=lambda) %>% autoplot

autoplot(cangas)

retaildata <- read.csv("retail.csv")
mytimeseries <- ts(retaildata[,4], frequency=12, start=c(1982,4))
(lambda <- BoxCox.lambda(mytimeseries))

## [1] -0.3912717

mytimeseries %>% BoxCox(lambda=lambda) %>% autoplot()

Lab Session 11

bricksq %>% mstl(s.window=11) %>% autoplot()

bricksq %>% mstl(s.window=11) %>% seasadj() %>% autoplot()

bricksq %>% mstl(s.window=11, robust=TRUE) %>% seasadj() %>% autoplot()

Lab Session 12

sa <- bricksq %>% mstl(s.window=11) %>% seasadj()
sa %>% naive() %>% autoplot()

fc <- stlf(bricksq, s.window=11, method='naive')
autoplot(fc)

checkresiduals(fc)

Lab Session 13

fit <- mstl(mytimeseries, lambda=lambda)
autoplot(fit)

mytimeseries %>% stlf(lambda=lambda) %>% autoplot()

Now try it on the training/test split

train <- window(mytimeseries, end=c(2010,12))
test <- window(mytimeseries, start=2011)

fc <- stlf(train, lambda=lambda)
autoplot(train) +
  autolayer(fc) +autolayer(test, series="Tests")

accuracy(fc, test)

##                      ME      RMSE       MAE        MPE     MAPE      MASE
## Training set -0.3506253  8.334542  5.215022 -0.1703858 2.904052 0.2696797
## Test set     -9.7437987 13.250443 10.684962 -4.9152571 5.367077 0.5525417
##                    ACF1 Theil's U
## Training set 0.03921138        NA
## Test set     0.23659054 0.6344392

Lab Session 14

usnetelec %>% autoplot

usnetelec %>% diff %>% autoplot

usgdp %>% autoplot

usgdp %>% diff %>% autoplot

lambda <- BoxCox.lambda(mcopper)
mcopper %>% BoxCox(lambda=lambda) %>% autoplot

mcopper %>% BoxCox(lambda=lambda) %>% diff(lag=12) %>% autoplot

enplanements %>% log %>% autoplot

enplanements %>% log %>% diff(lag=12) %>% autoplot

enplanements %>% log %>% diff(lag=12) %>% diff %>% autoplot

visitors %>% autoplot

lambda <- BoxCox.lambda(visitors)
visitors %>% BoxCox(lambda=lambda) %>% autoplot

visitors %>% BoxCox(lambda=lambda) %>% diff(lag=12) %>% autoplot

visitors %>% BoxCox(lambda=lambda) %>% diff(lag=12) %>% diff %>% autoplot

mytimeseries %>% 
  BoxCox(lambda=lambda) %>%
  diff(lag=12) %>%
  diff(lag=1) %>%
  autoplot

Lab Session 15

wmurders %>% autoplot()

wmurders %>% log() %>% autoplot()

fit <- auto.arima(wmurders, lambda=0, stepwise=FALSE, approximation=FALSE, max.p=10, max.q=10, max.order=10)
checkresiduals(fit)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,2,1)
## Q* = 10.562, df = 8, p-value = 0.2278
## 
## Model df: 2.   Total lags used: 10

forecast(fit) %>% autoplot()

wmurders %>% ets() %>% forecast() %>% autoplot()

wmurders %>% ets(lambda=0, model="AAN") %>% forecast() %>% autoplot()

Lab Session 16

lambda <- BoxCox.lambda(train)
auto.arima(train, lambda=lambda)

## Series: train 
## ARIMA(1,1,0)(1,1,0)[12] 
## Box Cox transformation: lambda= -0.2741334 
## 
## Coefficients:
##           ar1     sar1
##       -0.2188  -0.3834
## s.e.   0.0539   0.0510
## 
## sigma^2 estimated as 0.0001804:  log likelihood=960.12
## AIC=-1914.23   AICc=-1914.16   BIC=-1902.82

(arimamod <- auto.arima(train,
              D=1, lambda=lambda,
              stepwise=FALSE,
              approximation=FALSE))

## Series: train 
## ARIMA(3,1,0)(2,1,0)[12] 
## Box Cox transformation: lambda= -0.2741334 
## 
## Coefficients:
##           ar1      ar2     ar3     sar1     sar2
##       -0.2192  -0.0053  0.1307  -0.5164  -0.3482
## s.e.   0.0552   0.0559  0.0545   0.0527   0.0523
## 
## sigma^2 estimated as 0.0001574:  log likelihood=982.76
## AIC=-1953.51   AICc=-1953.26   BIC=-1930.68

(etsmod <- ets(train))

## ETS(M,Ad,M) 
## 
## Call:
##  ets(y = train) 
## 
##   Smoothing parameters:
##     alpha = 0.8031 
##     beta  = 0.0049 
##     gamma = 1e-04 
##     phi   = 0.9767 
## 
##   Initial states:
##     l = 63.945 
##     b = 0.2145 
##     s = 0.9895 0.9373 1.0084 1.1736 1.0097 1.0174
##            0.9782 0.9979 0.9931 0.9466 0.9759 0.9725
## 
##   sigma:  0.0457
## 
##      AIC     AICc      BIC 
## 3383.427 3385.526 3452.611

f1 <- snaive(train, h=length(test))
f2 <- hw(train, h=length(test), seasonal='multi')
f3 <- forecast(etsmod, h=length(test))
f4 <- stlf(train, lambda=lambda, h=length(test))
f5 <- forecast(arimamod, h=length(test))

checkresiduals(f1)

## 
##  Ljung-Box test
## 
## data:  Residuals from Seasonal naive method
## Q* = 1256.8, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

checkresiduals(f2)

## 
##  Ljung-Box test
## 
## data:  Residuals from Holt-Winters' multiplicative method
## Q* = 43.167, df = 8, p-value = 8.171e-07
## 
## Model df: 16.   Total lags used: 24

checkresiduals(f3)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(M,Ad,M)
## Q* = 40.992, df = 7, p-value = 8.126e-07
## 
## Model df: 17.   Total lags used: 24

checkresiduals(f4)

## 
##  Ljung-Box test
## 
## data:  Residuals from STL +  ETS(M,A,N)
## Q* = 91.221, df = 20, p-value = 4.531e-11
## 
## Model df: 4.   Total lags used: 24

checkresiduals(f5)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(3,1,0)(2,1,0)[12]
## Q* = 54.96, df = 19, p-value = 2.357e-05
## 
## Model df: 5.   Total lags used: 24

accuracy(f1,test)["Test set","RMSE"]

## [1] 23.40458

accuracy(f2,test)["Test set","RMSE"]

## [1] 59.85277

accuracy(f3,test)["Test set","RMSE"]

## [1] 25.43061

accuracy(f4,test)["Test set","RMSE"]

## [1] 34.82789

accuracy(f5,test)["Test set","RMSE"]

## [1] 105.5313

autoplot(f5) +
  autolayer(test, series="New data")

Lab Session 17

(etsmodel <- ets(mytimeseries))

## ETS(M,Ad,M) 
## 
## Call:
##  ets(y = mytimeseries) 
## 
##   Smoothing parameters:
##     alpha = 0.7718 
##     beta  = 0.0025 
##     gamma = 0.0561 
##     phi   = 0.98 
## 
##   Initial states:
##     l = 64.0899 
##     b = 0.2813 
##     s = 0.9869 0.9362 1.0043 1.1705 1.0216 1.0293
##            0.9848 0.9998 0.993 0.9436 0.9694 0.9606
## 
##   sigma:  0.0464
## 
##      AIC     AICc      BIC 
## 4404.041 4405.697 4477.272

accuracy(etsmodel)

##                     ME     RMSE      MAE       MPE     MAPE     MASE
## Training set 0.5339249 9.613425 6.241919 0.1970991 3.257992 0.333491
##                     ACF1
## Training set -0.03180379

etsfc <- function(y,h)
{
  y %>% ets(model="MAM", damped=FALSE) %>% forecast(h=h)
}
e <- tsCV(mytimeseries, etsfc)
sqrt(mean(e^2, na.rm=TRUE))

## [1] 10.62937

ggtsdisplay(e)

## Warning: Removed 1 rows containing missing values (geom_point).

lambda <- BoxCox.lambda(mytimeseries)
(arimamodel <- auto.arima(mytimeseries, lambda=lambda))

## Series: mytimeseries 
## ARIMA(0,1,1)(2,1,1)[12] 
## Box Cox transformation: lambda= -0.3912717 
## 
## Coefficients:
##           ma1    sar1     sar2     sma1
##       -0.2359  0.0307  -0.1728  -0.8660
## s.e.   0.0456  0.0615   0.0581   0.0476
## 
## sigma^2 estimated as 3.496e-05:  log likelihood=1548.16
## AIC=-3086.32   AICc=-3086.17   BIC=-3066.13

accuracy(arimamodel)

##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -0.3994054 9.337973 6.191385 -0.2677669 3.232095 0.3307911
##                     ACF1
## Training set 0.004237054

r lab17d} arimafc <- function(y,h) { y %>% Arima(order=???, seasonal=???, lambda=??) %>% forecast(h=h) } e <- tsCV(mytimeseries, arimafc) sqrt(mean(e^2, na.rm=TRUE)) ggtsdisplay(e)

etsfc <- function(y,h)
{
  y %>% ets(y) %>% forecast(h=h)
}
arimafc <- function(y,h)
{
  y %>% auto.arima(y, lambda=BoxCox.lambda(y)) %>% forecast(h=h)
}

e1 <- tsCV(mytimeseries, etsfc)
e2 <- tsCV(mytimeseries, arimafc)
#e <- cbind(e1, e2) %>% na.omit()
#colMeans(e^2)